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Abstract. Little is known about morphological instability of a solidification front 
during the crystal growth of a thin film of fiowing supercooled liquid with a free surface: 
for example, the ring-like ripples on the surface of icicles. The length scale of the ripples 
is nearly 1 cm. Two theoretical models for the ripple formation mechanism have been 
proposed. However, these models lead to quite different results because of differences 
in the boundary conditions at the solid-liquid interface and liquid-air surface. The 
validity of the assumption used in the two models is numerically investigated and 
some of the theoretical predictions are compared with experiments. 
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1. Introduction 

A thin liquid film flowing down a rigid wall is often observed in everyday life. A large 
number of studies on the instability of a viscous liquid layer running down a wall 
have been done (Benjamin 1957, Oron et al 1997). However, little is known about 
the morphological instability of the solid-liquid interface during the crystal growth of 
a flowing liquid film: for example, the ring-like ripples on icicles as shown in figure [1] 
(a). Icicles grow when their surfaces are covered with a supercooled water film and 
the latent heat is released through the water film into the ambient air below °C 
(Makkonen 1988). It is well known that the ripples on the surface of icicles have a very 
regular spacing of about 9 mm (Maeno et al 1994). A pattern similar to ripples on 
icicles was experimentally produced on the surface of a wooden round stick and that of 
a gutter on an inclined plane, by supplying water from their top in a cold room below 
°C (Matsuda 1997). He found that these also have centimeter-scale ripples on their 
surfaces. 

What determines the length scale of ripples on the surface of icicles? Although this 
has been a familiar phenomenon for people in cold regions for a very long time (Terada 
1947), nobody has been able to explain the details until recently. A first theoretical 
attempt to explain ripple formation on icicles was made in (Ogawa and Furukawa 2002). 
After that a quite different ripple formation mechanism was proposed by one of the 
authors (Ueno 2003, Ueno 2004, Ueno 2007). Although the governing equations in both 
papers were basically identical, completely different results were obtained due to some 
differences in boundary conditions. The author (Ueno 2003, Ueno 2007) predicted that 

(i) the wavelength of ripples increases with a decrease in the angle of the inclined plane, 

(ii) the wavelength increases only gradually with an increase in the water supply rate per 
width, and (iii) the ripples move upward at about half speed of the mean growth rate 
of icicle radius. The experimental results by Matsuda were obtained at a fixed water 
supply rate. We conducted similar experiments for various inclination angles and water 
supply rates. The purpose of this paper is to numerically investigate the validity of the 
assumptions used in both models (Ogawa and Furukawa 2002) and (Ueno 2003), and to 
compare the above theoretical predictions (i), (ii) and (iii) with our own experimental 
results. 

2. Theoretical framework 

We consider an ice growth on a substrate by supplying water from the top, as shown 
in figure [1] (b). One side of the water film is a water-air surface and the other side is 
growing ice. As a result of the instability of the ice-water interface as shown in figure [1] 
(c), the fiow in the water film can be changed depending on the morphology of ice. 
In previous papers (Ueno 2003, Ueno 2004, Ueno 2007), we assumed a semi-infinite 
ice layer and no airfiow ahead of the water-air surface. In this paper, we extend the 
previous theoretical framework to include heat conduction from the ice-water interface 
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Figure 1. (a) Ripples on natural icicles hanging from a roof (Feb 2009, Chicoutimi, 
Canada), (b) and (c) are schematic views of ice growth from a thin film of supercooled 
water flowing down a substrate inclined at angle to the horizontal. Flow is driven 
by gravity, and g is the gravitational acceleration, (b) shows an unperturbed state 
of the ice-water interface and water-air surface, Kq and uq are the thickness and the 
surface velocity of the flowing water fllm. Tgubi Tgi, Tia and Too are temperatures at 
the substrate {y = — ^o), ice- water interface (y = 0) and water-air surface {y = ho) and 
y = ho + 6o, respectively. The temperatures of the ice Tg, supercooled water Ti and air 
T(j are below °C. V" is an unperturbed ice growth rate, (c) shows a perturbed state 
of the ice-water interface and water-air surface. 

into the substrate thorough the ice with finite thickness. Not only does this theoretical 
framework enable us to rewrite the governing equations and boundary conditions in a 
more tractable form to solve numerically, but it also let us easily compare the difference 
between the models (Ogawa and Furukawa 2002) and (Ueno 2003). 

Instead of dealing with the complete geometry of the icicle, round stick and gutter 
on an inclined plane, the theoretical analysis is assumed to be restricted to two- 
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Table 1. Symbols. 


symbol 


aennition 


Q/l 


water supply rate per width 


e 


angle of inclined plane with respect to horizontal 


bo 


thickness of growing ice 


ho 


thickness of unperturbed water film 


5o 


characteristic length scale of thickness of thermal boundary layer in air 


Uo 


surface velocity of flowing water film 


c 


disturbed ice-water interface 




disturbed water-air surface 


-^sub 


temperature of substrate 


Tsi 


temperature at ice-water interface 


Tla 


temperature at water-air surface 


T 

oo 


ambient air temperature 


AT,; 


temperature deviation from T,/ 


AT,, 


temperature deviation from Tia 


Rei 


Reynolds number of water film flow 


Pei 


Peclet number of water film flow 


fi 


non-dimensional amplitude of perturbed stream function in water film 


Hi 


non-dimensional amplitude of perturbed temperature in water film 




non-dimensional wave number 


a 


restoring force due to surface tension and gravity 


(r) 


non-dimensional amplification rate of disturbed ice-water interface 




non-dimensional translational velocity of disturbed ice-water interface 




ratio of thermal conductivity of ice to that of water 




ratio of unperturbed temperature gradient in ice to that in water 




phase shift of water-air surface against ice-water interface 




phase shift of temperature at ice-water interface against ice-water interface 




phase shift of temperature at water-air surface against ice-water interface 




phase shift of heat flux at ice-water interface against ice-water interface 



dimensional vertical cross-sections of their objects, as shown in figures [T] (b) and (c). 
The X axis is parallel to semi-parabolic shear flow direction, and the y and y axes are 
normal to it. y is a laboratory frame, and y is a moving frame with an undisturbed ice 
growth rate V. and Uq are the thickness and the surface velocity of an undisturbed 
flowing supercooled water film. For typical values of water supply rate per width Q/l 
(Maeno et al 1994), mq ~ 1 cm/s and ho ~ 100 /im. Since actual thickness of the water 
film is very thin, the thickness of water film is drawn exaggeratedly in figures [1] (b) and 
(c). For convenience, we list only non-standard or particularly important symbols in 
Table [H 
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2.1. Governing equations 

The velocity components ui and vi in the x and y directions in the water film flowing 
down an inclined plane at angle 6 with respect to the horizontal are governed by the 
Navier-Stokes equations driven by gravity and the continuity equation (Landau and 
Lifshitz 1959): 

dui dui dui I dpi ( d'^ui d'^ui\ . 

^ + ^'^ + ^'^ = + ^i\^rT + ^TT ] + 9sm9, (1) 

ot ox oy pi ox ' '^""^ '^"■^ ' 

dvi dvi dvi 1 dpi 
dt dx dy pi dy 

dui ^ dvi 
dx dy 




dvi dvi dvi I dpi , ^ ^ , 

^ + ^'^ + ^'— = -T— + ^^^1^ + 1^:^) (2) 

dui dvi „ 

^ + = 0, (3) 



where t is time, pi the pressure, pi = 1.0 x 10'^ kg/m^, the density of water, ui = 1.8 x 10^^ 
m^/s, the kinematic viscosity of water, g=9.8 m/s^, the gravitational acceleration. From 
([3]), using the stream function ipi, ui and vi can be expressed as ui = dipi/dy and 
vi = —dipi/dx. 

The equations for the temperatures in the ice Tg, water T/, and air Ta are (Landau 
and Lifschitz 1959, Caroli et al 1992) 

^^'^ = J + ) , (6) 




where Kg = 1.15 x 10~^ m^/s, = 1.33 x 10"'' m^/s and = 1.87 x 10"^ m^/s are the 
thermal diffusivities of the ice, water and air, respectively. We can neglect the second 
terms on the left hand side of (jl]) and ([6]) because ice grows very slowly (Ueno 2003). 

2.2. Boundary conditions 

2.2.1. Hydrodynamic boundary conditions Neglecting the density difference between 
ice and water, both velocity components at a disturbed ice-water interface must satisfy 
(Myers et al 2002a, Ogawa and Furukawa 2002) 

Ul\y=C = 0, Vl\y=^ = 0. (7) 

Except for ([7]), the following boundary conditions are the same as those used in the 
stability analysis of a viscous liquid layer flowing down a rigid wall (Benjamin 1957). 
The kinematic condition at a disturbed water-air surface is 
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At the water-air surface the shear stress must vanish: 
dui 
dy 



dvi 
H 

y=i dx 



0, 



and the normal stress including the stress induced by the surface tension 7 
N/m of the water-air surface must balance the atmospheric pressure Pq: 

-3/2 



I dvi 

Pl\y=^ + 2pil^l — 
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-Pn. 



(9) 

7.6 X 10-2 
(10) 



2.2.2. Thermodynamic boundary conditions In the model (Ogawa and Furukawa 
2002), the continuity of the temperature at a disturbed ice- water interface, y = ({t,x), 



IS 



Ts\y=(: - Ti\y=i^ — Tsl. (11) 

If we can neglect the Gibbs-Thomson effect (temperature depression due to the curvature 
of the solid-liquid interface) (Caroli et al 1992), T,/ is assumed to be the equilibrium 
freezing temperature (Tgi = °C for pure water). On the other hand, in the model 
(Ueno 2003), the continuity condition is represented as follows: 

T,\y=^^Ti\y=^^T,i + ^T,i. (12) 

We will discuss in Section 4 that the temperature at a disturbed ice-water interface 
under a thin shear flow is not necessarily a constant T^;, but that it deviates by AT^; 
from Tsi- The heat conservation at the ice- water interface is (Langer 1980, Caroli et al 
1992, Ueno 2003) 




K 



dy 



y=C dy 



y=C 



(13) 



where L = 3.3 x 10*^ J/m^ is the latent heat per unit volume, and Kg = 2.22 J/(mKs) 
and Ki = 0.56 J/(m K s) are the thermal conductivities of the ice and water, respectively. 

In the model (Ueno 2003), the continuity of the temperature at a disturbed water- 
air surface, y = ^{t,x), is 

Tl\y=^ ^ Ta\y=^ ^ Tia. (14) 

We will discuss in Section 4 that the temperature at a disturbed water-air surface should 
remain at a constant Tia, which will be determined from the continuity of heat flux at 
the water-air surface. On the other hand, in the model (Ogawa and Furukawa 2002) the 
continuity of the temperature is 



Ti\y=(^ — Ta\y=(^ — Tia + ATia, 



(15) 



which means that the temperature at a disturbed water-air surface deviates by ATia from 
Tia- The heat conservation at the water-air surface is given by (Ogawa and Furukawa 
2002, Ueno 2003) 

dTi 

— K 



-K, 



dy 



y=i " dy 



(16) 
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where Ka = 0.024 J/(mKs) is the thermal conductivity of the air. 

We will see later that completely different results in the two models arise from the 
different boundary conditions between ( fTTl) . ( fT5i) in (Ogawa and Furukawa 2002) and 
(fT2|) . (Ill]) in (Ueno 2003). AT^; in (fT2|) and AT;^ in (ITSl) cannot be determined a priori, 
but will be determined after solving the equation for the temperature in the flowing 
water film. 



2. 3. Equations and solutions for unperturbed and perturbed fields 



Since a ring-like structure encircles the icicles and there is no noticeable azimuthal 
variation on the surface of the icicles (see figure [1] (a)), it is sufficient to consider 
only a one dimensional perturbation in the x direction of the ice-water interface, 
C{t, x) = exp[at+ikx], where k is the wave number and a = a^'^''+ia^^\ with a*^''^ being 
the amplification rate and Vp = —a^^'^/k being the phase velocity of the perturbation, 
and (k is a small amplitude of the ice-water interface. We separate ^, iph Ph Ts, Ti and 
Ta into unperturbed steady fields and perturbed fields with prime as follows: ^ = /io + ^', 
tpi = ipi + ip'i, pi = Pi + p'l, Ts = fs + T^, Ti = fi + T/ and Ta = fa + T^. We suppose 
that the respective perturbed parts are expressed as follows: 



1 ^'it,x) \ 




( ) 






Fi{y) 


p'i{t,x,y) 




My) 


Ts{t,x,y) 




9s{y) 


Tl{t,x,y) 




9i{y) 


\ T'ait,x,y) ) 




\ 9a{y) J 



exp[crt + ikx], 



[17) 



where ^k, ^i, ^i, ds, Qi and Qa are the amplitudes of respective perturbations and 
they are assumed to be of the order of Cfc- The calculation in the previous paper 
(Ueno 2003) was based on a linear stability analysis taking into account only the 
first order of Cfc- Furthermore, two approximations were used. The first is the long 
wavelength approximation (Benjamin 1957, Oron et al 1997), which is valid when the 
water film thickness is much less than the characteristic length scale of ripples. Defining 
a dimensionless wave number by = kho, we neglected the higher order of /i. The 
second is the quasi-steady state approximation (Langer 1980, Caroli et al 1992). We 
neglected the time derivative term of m;, vu Tg, Ti, Ta, ^ in ([1]), ([2]), (j4]), ([5]), IQ and 
(IE]) because these fields respond relatively rapidly to slow development of the ice-water 
interface perturbation. In order to check numerically the validity of the analytical results 
obtained under the long wavelength approximation in the previous papers (Ogawa and 
Furukawa 2002, Ueno 2003), we retain the higher order of fi in the following perturbed 
parts. 

The equations of the unperturbed part in ([T]) and ([2]) are, respectively, 

0, -^-gcos9 = 0. (18) 



dy2 



Pi 
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With the no-shp condition at the ice-water interface, Ui\y=o = 0, the free shear stress at 
the water-air surface, dUi/dy\y=ho = 0, and Pi\y=ho = Pq, the solutions are 

Ui{y) = uo |2|^ - (l^)'l , Pi{y) = Po- PigcosOiy - ho), (19) 

where uq = K^gsuiO / {2vi) is the surface velocity of the water film. In the absence of ice 
growth, the water supply rate per width is given by Q/l = Jq'^ Ui{y)dy = 2Mo/io/3 in an 
undisturbed state (Benjamin 1957, Landau and Lifschitz 1959), from which ho and Uq 
can be expressed with respect to experimentally controllable parameters Q/l and 6 as 
follows: ho = [3i/,/(^sin0)]V3(Q//)i/3 and % = [9g sme/{8ui)Y/\Q/lf/\ 

From the dimensional consideration and the assumption that Fi in f lT7|) is of 
the order of (k, we assume Fi{y) = uofi{y)(k- Substituting Fi and 11/ in ( |T7|1 
into the perturbed part of ([1]) and ([2]), and finally eliminating 11; from them by 
cross differentiation, we obtain the Orr-Sommerfeld equation for the non-dimensional 
amplitude of the perturbed stream function fi (Benjamin 1959): 

where = y/ho, = Ui/uq = 2y^ — y1 and Rei = uoho/ui = 3Q/{2lh'i) is the Reynolds 
number. 

Linearizing the boundary conditions (17|)- ffTU]) at the unperturbed ice- water interface 
y = and water air surface y = ho, the perturbed part of the equations are, respectively, 
dUi/dy\y=oC + u\\y=o = 0, vl\y=o = 0, Ui\y=hod^' / dx = v'i\y=ho, (i'^Ui/dy'^\y=ho^' + 
du'i/dy\y=ho + dv'i/dx\y=ho = and -{dPi/dy\y=h^^' + p'i\y=ho) + 2piPidv[/dy\y=ho - 
'-yd'^^' /dx"^ = 0. Using u'l = dip'Jdy and v[ = —dipl/dx, the above equations can be 
expressed as respectively (Ueno 2003): 





fl\y^=0 — 0, fl\y, = l(k — —Ul:t\y^ = l^k, 



dm. 



- (i/i/?e/f7i*|j/,=i + 3/i^)^ +ifiRei'^^ Ji\y,=]\ Ck = '^aik- (21) 
dy^ s/.=i dy^ y*=i dy^ y*=^ J 

While deriving the last equation in (j2Tl) . we have used 11/ = piu^ / ho{l / {iiiRei){d? fi / dyl— 
ji^dfi/dy^) — Ui^dfi/dy^ + dUi^/dy^fi}C,k obtained from the perturbed part of pressure 
gradient term in ([1]). Here 

2 / - ^ 2 



a = 2(cote)/x + -^(-l) /x^ (22) 

in the last equation in (1211) is the parameter to characterize the effect of surface tension 
and gravity on the water-air surface, which was referred to as the restoring force in 
the papers (Benjamin 1957, Ueno 2003). a = [j/ipig)Y^^ is the capillary length 
associated with the surface tension 7 of the water-air surface (Landau and Lifschitz 
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1959). The third equation in f l2ip gives the relation between the amphtude of the ice- 
water interface, (k, and that of the water-air surface, C,k- In the case of a disturbed 
ice-water interface and water-air surface too, ui{y)dy = J^{Ui{y) + u'i{x,y)}dy = 
2uoho/3 + Uo{^k + /«|i/*=iCa:) exp[(Tt + ikx] up to the first order of (k must be equal to 
Q/l, from which we again obtain the relation ^k = —fi\y,=iCk- Noting that Uit,\y^=i = 1, 
dtJiJ dy^\y^=o = 2 and d'^UiJ dyl\y,^=i 
four boundary conditions to solve ([20 

dfi 



-2 and using C.k = -fi\y,=iCk, (EI]) leads to 



+ 2 = 0, 



fi 



dyl 



dyl 



- {ifiRei + 3/i )- — 
y.=i dy^ 



= 0. (23) 
yields 



Linearizing ( [T2|) and (fT3|) at the unperturbed ice-water interface y 
respectively, to the zeroth order in (k, 



T, 



l\y=0 



sh 



dbo 



and to the first order in (k, 
"dT, 



AT,, 



dy 



y=0 



Ck + g. 



s\y=0 



LaQk = Ks — 
dy 



exp[c7t + ikx] 



df, 



df. 



dy 



y=0 



K, 



d7] 
dy 



y=o' 



(24) 



dy 



y=0 



Ck + 9i\ 



y=0 



exp[crt + ikx], 



y=0 



y=o " dy 

Next, linearizing f|T5|) and f|T6|) at the unperturbed water-air surface y 
respectively, to the zeroth order in ^k, 



T, 



l\y=ho 



a\y=ho 



T, 



lai 



-K, 



and to the first order in C,k, 
'dfn 



AT, 



la 



dy 



y=ho 
Ki 



ik + gi\ 



y=ho 



'd'T, 



dy' 



y=ho ay 



exp[crt + ikx] 
dgi\ ^ 



dy 
df„ 



y=ho 



-K, 



dfa 



dy 



y=ho 



y=ho , 



d^T. 



dy' 



dy 

+ ga\y=ho 

. dga 
y=ho ay 



y=ho 



(25) 
ho yields 

(26) 



exp [at + ikx] , 



y=ho 



(27) 



Substituting Tg = Tg + gsexp[at + ikx], = + (7/ exp [at + ikx] and = 
Ta + (^a exp [at + ikx] into (jl]), ([5]) and (Q, the equations for Tg, Ti, Ta and gs, gi, 
ga are obtained. With the following boundary conditions: Ts\y=^bo = ^sub, ^s|y=o = 
fi\y=o = Tsi, fi\y=ho = Ta\y=ho = Tia and fa\y=ho+5o = Too, as shown in figure [D (b), we 
obtain linear temperature profiles 

fsiy) = Ts^b + Gsiy + bo), fi{y) = T^i - G,y, f„(y) = T,, - Ga{y - h), (28) 

where Gs = (Tsi - Tsub)/&o, Gi = {Tsi - Tia)/ho and Ga = {Tia - Too)/5o are the 
unperturbed part of temperature gradient. Here we have assumed that air temperature 
is approximately Too at a length scale 5q. In the presence of airflow, heat transport is 
greatly influenced by the convection, and is then regarded as a characteristic length 
scale of the thickness of the thermal boundary layer (Short et al 2006). 
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When we consider the icicle as a cyhndrical object, the solutions of the unperturbed 
velocity and temperature profiles in the water film in the cylindrical coordinate under 
the assumption of axial symmetry with the same boundary conditions as the planar case 
are [/;,(r,) = Rl{l + l/ R:)Hn{rJ R;) -{rl- Rl)/2 andf,,(r,) = {Ti-T,{) / {T,i-Ti,) = 
—R^\Yi{r^./R^.), where r^, = r/ho and i?* = R/ho, r and R being the radial coordinate 
and the icicle radius, respectively. When we express Uu{r^) and T/*(r^,) with respect 
to using the relation = R^ + y^, the planar velocity and temperature profiles 
Ui*{y*) = 2y* — yl and Ti^.{y^) = —y^ are retrieved because y^/ R^, ^ 1 in the water film 
(0 < < 1) when the icicle radius R is much greater than the thickness of water film 
/lo, i.e., i?* ^ 1. That is why the icicle geometry was approximated in the Cartesian 
coordinates. 

Substituting the solutions of % and Ta in ( l28l) into the second equation in ( l26l) . Tia 
in and ( |T5l) is obtained as follows: 

rp _|_ Kg hp rp 

Next substituting the solutions of Tl, and Ti in (128|) into the second equation in (12^ and 
using (j2H]) yields an unperturbed ice growth rate approximately: 

_ d6o _ Ks Tsi — Tsub Ka Tsi — Tqq 
^=^~T bo +X 5o ' ^ ^ 

because Ka/Ki <^ 1 and /iq/^o ^ 1. In the presence of airflow, (1291) and (1301) are 
replaced by T,, = {%i + (ir,/irO(/^o/5o)G',*Too}{l + (ira/irz)(/io/5o)G',4 and dh/dt = 
{Ks/ L){Tsi-Tsuh) /bo + {Ka/ L){Tsi-T^)GaJ5o, where Ga* = -Sq/ {Tia-T^)dTa/dy\y=ho 
is the dimensionless air temperature gradient at the unperturbed water-air surface and 
depends on the Prandtle number of the air. In order to estimate the value of Ga*, we 
need the exact temperature distribution by solving the coupled Navier-Stokes and 
heat transport equations in the air (Short et al 2006). As we will see later, however, as 
the air temperature gradient Ga at the water-air surface does not affect the wavelength 
of ripples on icicles, we assumed a linear air temperature proflle ahead of the water-air 
surface, i.e., we consider the case of Ga* = 1- 

Figure [2] shows the ice thickness determined by integrating (1301) . subject to 6o = at 
time t = 0, for different temperatures Tg^b of a substrate. If heat conduction through the 
ice is negligible, is proportional to the time t and V = 2.6 mm/h for the parameters 
shown in figure [2l This value is of the same order as our experimental measurement of 
ice growth rate (see Section 5). When T^uh < °C, is proportional to t^/^ because the 
first term on the right hand side in (|30ll is dominant while the thickness of ice is small. 

From the dimensional consideration from the first equation in fl25l) . we assume 
gi{y) = Hi{y)Gi(k- Then the perturbed part of ([5]) yields the equation for the non- 
dimensional amplitude Hi (Ueno 2003): 

^ = (/x^ + i^PeiUu)Hi - i/^Pei^fi, (31) 
dy: dy* 
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V 1 4 £ I in 



Figure 2. Change in ice thickness with time for different temperatures Tsub of a 
substrate. Sq = 1.0 mm (Short et al 2006), T,/ = °C and -10 °C. 



where Ti^{y*) = {Ti{y^) — Tsi)/{Tsi — Tia) = —y* is the dimensionless temperature profile 
of the water film in the unperturbed state, and where Pei = Uoho/Ki = 3Q/{21k,i) is the 
Peclet number defined as the ratio of the heat transfer due to the water flow to that 
due to the thermal diffusion in the water film. 

In the case of ATia = as in (fT^ (Ueno 2003), the first equation in (1271) gives 
Hi\y=hoCk = and ga\y=ho = Ga^k- Then the solution of the equation d^ga/dy"^ = k'^Qa 
with the boundary conditions of ga\y=ho = Gaik and ga\y=oa = is given by 

gaiy) = exp[-k{y - hQ)]Ga^k- (32) 

Substituting ( 132|) into the second equation in ( 1271) and using the second equation 
KiGi = KaGa in ( 126|) . yields dHi/dy\y=ho(k = —k^k- Accordingly, using the relation 
^k = —fi\yt=iCk the boundary conditions to solve ( pTl) are 

m 

dy* 

Using the solution Hi satisfying the boundary conditions in (133|) . from the first equation 
in ( 125|) . ATsi and gs\y=o are finally determined to the first order in (j.: AT^i = 
{Hi\y^=o-l)GiCkGwWt + ikx] aiidgs\y=o = i-Gs/Gi + Hi\y=o-l)GiCk,T^Gspective\Y. The 
solution to the equation d'^gs/dy'^ = k'^gs with the boundary conditions of gs\y=~bo = 
and gs\y=o = {-Gs/Gi + Hi\y=o - l)Gi(k is 



Hl\y^ = l — -fl\y, = l, --J— ^ , — -p-fl\y, = l- (33) 



\ Gi J smh(fc6o) 

On the other hand, in the case of ATsi = as in ( ITTl) (Ogawa and Furukawa 2002), 
the first equation in ( 125|) gives gs\y=o = —GsCk and Hi\y=Q = 1. Then the solution 
to the equation d'^gs/dy'^ = k'^gs with the boundary conditions of gs\y=-bo = and 
gs\y=Q = -GsCk is given by 

smh[k{y + bp)] - 

9^iy) = ^TTTin — 35 

smh(A;oo) 
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Illy* 
HI ^ 



[fl\y,=, - ^my*=l + fl\y,=l)] ■ (37) 



Furthermore, the first equation in 027]) gives ga\y=ho = GaCk + Gi{Hi\y=hoCk - Ck)- 
Then the solution to the equation d'^Qa/dy'^ = k'^Qa with the boundary conditions of 
ga\y=ho = Gaik + G i{H i\y=hQC,k " ^k) and ga\y=oo = takes the form: 

9aiy) = exp[-A;(?/ - ho)]{Ga^k + Gi{Hi\y=hoCk - 6)}, (36) 

where Hi is not yet determined. Substituting (!36|) into the second equation in (1271) 
yields dHi/dy\y=hoKiGi(k = -kKa{Ga^k + Gi{Hi\y=hoCk-^k)}- Then, using the relation 
= —fi\y^=iCk the boundary conditions to solve fl3Tl) are 

= —II, < f/ L. _i — 

Since Ka/Ki <^ 1, we can neglect the second term of the second equation in (IHTIl 
and the result is the same as the second equation in fl33|l . Finally, substituting the 
solution of Hi satisfying the boundary conditions in fl57|) into fl5Bl) . the solution of Qa is 
determined and AT;^ is obtained from the first equation in (!27|) to the first order in C,k: 
AT/a = [Hi\y^=i + //|j^^=i)GzCfc exp[(Tt + i/cx]. We will see later that the difference in the 
boundary conditions between (l33l) (Ueno 2003) and (1371) (Ogawa and Furukawa 2002) 
is the main cause leading to different results. 

A small perturbation of the ice-water interface in the non-dimensional form can be 
rewritten the following way: 

= C* = (5blni[exp(crt + ikx)] = 6b(t) sin[A;(x — Vpt)], (38) 

where Sb = Ck/ho, 6bit) = 6bexp{a^^'H) , and Im denotes the imaginary part of its 
argument. From the amplitude relation = —fi\y^=iCk, the corresponding perturbation 
of the water-air surface with an infinitesimal amplitude 6t = ^k/ho is given by 

y^: = ^* = 1 + lm[5texp{at + ikx)] 

= 1 + \fi\y^=i\Sb{t) sm[k{x - Vpt) + e^J, (39) 

where = [{-ft''^\y,=i)' + {-fP\y,=i)V^ cosB^, = -//^\.=i/|/d,,=i| and 

sin 9^^ = —//*■' 6^, being a phase difference between the ice-water 
interface and the water-air surface. When /i is small, fi\y,=i = fi^^\y,=i + ifl^^\y^=i ~ 
-36/(36 + a^) + i(-6a)/(36 + a^) (see equation (73) in (Ueno 2003)). When //^^|y.=i 
acquires non-zero values, as shown in figure [3] (a), there exists a phase shift of the water- 
air surface against the ice-water interface. The parameter a in fl22|) depends on fi, hence 
the amplitude and phase of the water-air surface relative to the ice-water interface can 
change depending on the wave number of the ice- water interface (figures [3] (a) and[6](e)). 
Using the second equation in ( !33l) . the temperature gradient at the perturbed water-air 
surface can be expressed as —dTi/dy\y=^ = Gi{l — dHi/dy\y^=i() = Gi{l — fifi\y^=iC). 
Figured (b) shows the behaviour of the real and imaginary part of —^fi\y^=i with respect 
to /i. It is found that the perturbed part in —dTi/dy\y=(^ increases with an increase in /i, 
but that it is suppressed due to the restoring force acting on the water-air surface. This 
indicates that when the ice-water interface and water-air surface are coupled, even if the 
modes are sinusoidal in the linear stability analysis, the amplitude and phase shift of 
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Figure 3. For Q/l = 50 [(ml/h)/cm] and 9 = 7r/2, (a) shows the dependence of the 
non-dimensional amplitude of the perturbed steam function in the water film, |/;|y,=i|, 
and its real part — imaginary part on the non-dimensional wave 

number /i. (b) represents the behaviour of —^fi^\y,^i and —^fl^^\y,^i against ^. 
Here /i = 1.0 corresponds to the wavelength of 580 /im. 



the water-air surface against the ice-water interface significantly influence the perturbed 
part of temperature in the water film through the boundary conditions in ( 133|) . 

On the other hand, the effect of the restoring force on the water-air surface was not 
taken into account in the model (Ogawa and Furukawa 2002), i.e., fi\y,=i = —1, which 
means that the amplitude of the water-air surface is the same as that of the ice-water 
interface, and that there exists no phase shift of the water-air surface against the ice- 
water interface. Then, the perturbed part in —dTi/dy\y=^ = Gi{l + only increases 
with fi. This is also the main difference between the two models. 



2.4- Dispersion relation 

Substituting gi{y) = Hi{y)Gi(k and ( IMll into the second equation in ( l25ll and using 
KiGi = KaGa yields the dispersion relation for the perturbation of the ice-water 
interface: 

where Kf = Kg/Ki = 3.96 is the ratio of the thermal conductivity of ice to that 
of water and Gf = Gs/Gi = {Ki/ Ka){So/bo)(Tsi - Tsuh)/(Tsi - T^) is the ratio 
of the unperturbed temperature gradient at the ice-water interface in the ice to 
that in the water. The real and imaginary part of ( HOl) yield the non-dimensional 
amplification rate ai^'^ = a'^^^ / {KaGa/ (Lho)} and the non-dimensional phase velocity 



V, 



Vp/{KaGa/L) = -a^''^ /{kKaGa/L), respectively. 



(r) 



y*=o 



+ K^^{~Gt + Ht^\y,=o-l), (41) 
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^^ + ir,>ff«U=ol, (42) 



(r) (i) 

where Hi and if/ are the real and imaginary part of Hi. 

Equations (120|) and (I3T]) were solved analytically with the boundary conditions 
( l23l) and ( l33l) under the long wavelength approximation neglecting the higher order 
of fi, except for retaining the second term in a because of {a/hoY ^ 1 (Ueno 2003). 
By transferring the variable ?/* into z = 1 — y*, the general solution of ( l3Tl) can be 
expressed as Hi{z) = Ci4>i{z) + C24>2{z) + ifiPei Jq {(j)2{z)(f)i{z') - (f)i{z)(f)2{z')} fi{z')dz' , 
where Ci and C2 are unknown constants, (pi{z) and 02 (-2) are the homogeneous solutions 
of (ED): d^(f)/dz^ - {/i^ + i^PeiUu{z)}(t) = 0, where Uu{z) = 1- z'^. We note that in 
the z coordinate, z = and z = 1 are positions of the unperturbed water-air surface 
and ice-water interface, respectively. Neglecting the fi^ term in the above homogeneous 
equation, 4>i{z) and 02 (-2) can be expanded in terms of fj,Pei as follows (see APPENDIX 
in (Ueno 2003)): 

(f)i(z) = l + i( -z^ - —zA fiPei + (-—z^ + —z'' - — UPeif + ■ ■ ■ , (43) 
' V2 12 / V 24 360 672 J ^' ' ^ ' 



(pJz) = z + i[-z-^ z^ fiPei + z^ + z^ z-^ (l^PeiY + ■■■. (44) 

^ ^ ' V6 20 / V 120 2520 1440 J ^' ^ ' 

The boundary conditions in ( j33l) give Ci = —fi\z=o and C2 = —fifi\z=o, respectively, 
because 0i|2=o = 1; 02 1 2=0 = 0, d0i/d2;|2=o = and d02/d2;|2=o = 1- Thus we obtain 

Hi{z) = -fiU^o {0i(z) + fi<f)2{z)} + ifiPei r {02(^)01 (^') - 0i(^)02(z')} M^')dz'. (45) 

Jo 

Since 

5 239 
0i|,=i = 1 + i—ifiPei) - 7777— (yuPe;)^ + • • • , 



12'^ " 10080 

. — (uPei) , 

60^^ ^ 3360 



52U=i = 1 + i7^(/wPe/) - -^{fiPeiY 



dz 

d02 
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i-{fiPei) - —ifiPei)' + 

1. ^ . 17 



1 + i-(aPei) - ^^(aPeiY + ■ ■ ■ , (46) 
dz ^=1 4^^^ ^ 1440 u , K J 

the ratios of the second order term in fiPei to the first order one in 011^=1, 02 [2=1, 

d(pi/dz\z=i and d(l)2/dz\z=i are about 5.7 x 10~^/iPe;, 3.3 x 10~^/iPe;, 9.3 x 10~^/iPe/ 

and 4.7 X lO'^/iPe^, respectively. The second order terms in fiPei are not negligible 

when fiPei ~ 10. This is possible when Q/l ^ 300 [(ml/h)/cm] for the wavelength of 

ripples on icicles. However, the typical values of Q over icicles are on the order of tens 

of ml/h and their radii are usually in the range of 1 ~ 10 cm, the value of Q/l is in 

the range 10 ~ 100 [(ml/h) /cm] (Maeno et al 1994, Short et al 2006). As far as we are 

limited to such a range of Q/l, fiRei ^ 1 and /iPe; ~ 1 for the length scale of ripples on 

icicles. We can neglect the /iPe; term in ( 120|) and ( l23l) . This corresponds to neglecting 

the inertia term of the full Orr-Sommerfeld equation, then we can approximate (l20l) as 

follows: d'^fi/dyf = 0. The solution of this equation with the boundary conditions in 
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is given by fi{z) = (-6 + iaz + 6z'^ - iaz^)/{6 - ia) (Ueno 2003, Ueno 2007). As 
far as fiPei ~ 1, it is sufficient to consider up to the ffist order in fiPei in fH3l) and (jHj) 
because the second order terms in jj^Pei in ( H6l) are very small compared to the ffist 
order terms in fiPei as estimated above. 

Furthermore, if the substrate is not sufficiently cold or if there is not heat conduction 
through the substrate, the unperturbed part of heat conduction to the interior of 
the icicle is negligible (Makkonen 1988). The semi-infinite ice layer approximation of 
previous paper (Ueno 2003) corresponds to the case of 6o ^ 1- In these situations, 
G'l = 0, hence dH]) and (gS]) yield (Ueno 2003) 



cr. 



(r) 



/i |36 - |a(yuPe/)| - |a(/iPe,) _^ |36 - ^a(/iPe;)| - ^a(/iPei) - 



36 + ^2 



36 + a2 



(47) 



1 



/i {6a + 9{nPei)} - \Q^{^Pei) 



36 + a2 

jj |6a + y (/iPei)| + 6a - ^a^{^Pei) 



36 + a2 



0, al and v. 



p* 



are independent of the 



It should be noted that in the case of Gf 
unperturbed air temperature gradient Ga- 

On the other hand, Hi\y^=Q = Hi\z=i = 1 in the case of ATgi = and fi\y,=i = —I 
(Ogawa and Furukawa 2002) and noting that Ka/Ki <^ 1, the solution Hi with the 
boundary conditions in (137|) is given by if; (2;) = {(1 — iyuPe;J|2=i)0i(2;) +/i(0i|^=i02(-2) — 
02|.=i0i(^))}M|.=i + ifiPeiI{z), where I{z) = J,'{Mz)Mz') - M^)Mz')}fiWz' 
(Ueno 2004). This solution gives —dHi/dy^\y^=o = dHi/dz\z=i = /x/0i|^=i and in the 
0, gll) and (iSD yield 



case of Gf 



a 



(r) 



2 ' 



^fiPei 



jptf 



(49) 



(50) 



We notice that fH9l) and fl50l) are the same results as equations (77) and (81) in (Ogawa 
and Furukawa 2002), except /iPe; is defined as a in their paper. It is remarked that ( l49ll 
and ( l50l) are obtained by expanding 011^=1 up to the second order in /xPe^, in contrast 
to (HZD and (jiHl). 



3. Comparison of analytical with numerical results 

In spite of using many approximations in previous papers (Ogawa and Furukawa 
2002, Ueno 2003), the analytical calculations to solve the equations for /; and Hi 
with appropriate boundary conditions were very complex and cumbersome. Instead, 
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by decomposing fi and Hi into its real part // ,Hl , imaginary part fi ,Hi , we 
performed numerical studies to solve the ordinary differential equations for fi'^\ fi^\ 
hI^'' and hI^^ in equations (120|) and (13T]) . with boundary conditions (!23l) and (l33l) and 
without approximations as used to derive the analytical results mentioned above. 

First, in the case of 7^ 0, we have to consider heat conduction through a finite ice 
thickness into the substrate. Figure H] (a) shows the non-dimensional amplification rate 
cr* versus the non-dimensional wave number /i for different values of Gf at Q/l = 50 
[(ml/h)/cm] and 6 = 7t/2. In the range of < < 0.3, the wavelengths are longer 
than that in the case of = as the value of increases. We find that ai^^ < 
for all fi above Gf = 0.3. This means that if we choose the parameters in Gf = 
{Ki/Ka){So/bQ){Tsi—Ts^]j)/{Tsi—Too) to satisfy Gf > 0.3, ripples do not appear on the ice 
surface. Indeed, this is relevant as an experimental result that no ripples were observed 
on the ice grown on an planar aluminum substrate by supplying water from the top of 
the apparatus in a cold room at temperature below °C (Matsuda 1997). A smooth 
ice surface was produced on the aluminum substrate. Matsuda states that since the 
thermal conductivity of aluminum is about 100 times greater than that of ice, most of the 
latent heat released at the ice-water interface is conducted into the aluminum substrate 
through the ice. If there exists heat conduction through the substrate, the continuity 
of heat flux at the boundary between substrate and ice is KsuhiTsuh — T!3ubo)/4ub = 
KsiTsi — Tsuh)/bQ, where Ksuh is the thermal conductivity of a substrate, /sub is the 
thickness of the substrate, Tsuh and Tsubo are temperatures at the boundary between 
substrate and ice and that at the other side of the substrate, respectively. From this, 
^sub = {T^uho+iKs/ Ksnh){knh/ bo)Tsi} / {l+{Ks/ Ksnh) {huh /bo)} is obtained. Substituting 
this Tsub into the above Gf, we obtain Gf = {Ki/ Ka){6o/bo){Tsi-Tsubo) / (Ts;-Too)l/{1 + 
(i^'s/ii'sub)(/sub/^o)}- When the thickness of ice grown on the planar aluminum substrate 
satisfies the condition 69 ^ {Kg/ Ks^^,)lsnh, the above Gf can be approximated as 
Gf = {Ki/Ka){So/bo){Tsi — Tsuw)/{Tsi — Too). Moreover if the other side of surface 
of the aluminum substrate is exposed to ambient cold air, we assume Tgubo = Too. Then 
Gf = {Ki/Ka){So/bo) satisfies the condition Gf > 0.3 when 60 < 100 ^o- While the 
thickness of ice growing on the planar aluminum substrate by supplying water satisfies 
this condition, ripples would not appear on the ice surface. 

In the following discussions and the next sections, we will focus on the case of 
Gf = 0. The solid lines in figures H] (b) and (c) show ai^^ in ( HTl) and fp* in (H2ll 
against /i at Q/l = 50 [(ml/h)/cm] and 9 = 7r/2, respectively. The dashed and dashed- 
dotted lines in figures H] (b) and (c) represent the first and second term in ( HTl) and ( H2l) . 
respectively. In the case of ATgi = (Ogawa and Furukawa 2002), the second terms in 
( 14T1) and ( l42l) must be equal to zero because Hi^^\y_^=Q = 1 and -ff/*^|y,=o = from the 
first equation in ( 1371) . On the other hand, the contribution of ATsi 7^ in the model 
(Ueno 2003) appears in the second terms in fHTl) and fj48|) . Indeed, as shown by the 
dashed-dotted line in figures H] (b) and (c), the contribution of the second terms to the 
total values of ai^"* and Vp^, is not negligible. However, we can make an approximation in 
the case of figure H] (b) because the second term (dashed-dotted line) is smaller than the 
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Figure 4. Non-dimensional amplification rate ct; versus non-dimensional wave 
number fj, at Q/l = 50 [(ml/h)/cm] and 9 = 7t/2. (a) is for different values of Gf. 
(b) and (c) are in the case of Gf — 0. Solid lines: ai^'^ and Wp* are obtained from 
(|4T|) and (|42|). respectively; dashed lines: the contribution of the first term in (|4ll and 
(|42|) : daslied-dotted lines: the contribution of the second term in (|4T1) and p2)) . Here 
/i = 0.1 corresponds to the wavelength of 5.8 mm. 



first term (dashed line) and the wave number at which ai^^ acquires a maximum value is 
almost the same as that without the second term. When /i is small, we can approximate 
m as follows: ai'^ ^ dH^^ /dzU=, ^ -fifPU=o + {2/3)fiPeiffU=o-f^Pei Jo ff\z)dz. 
This indicates that the perturbed part of temperature gradient at the ice- water interface 
is affected by the amplitude //'^-'I^^q and //*^|^=o of the perturbed part of stream function 
at the water-air surface. Since /iPe^ ~ 1 and a ~ 1 for small /i and the typical range 
of Q/l = 10 ~ 100 [(ml/h)/cm] and 6 = tt/2, extracting the most dominant term from 
dHi^^ /dz\z=i and using fl22|) . we obtain 

^ 36,. - |a(/.PeO ^ ^ _ Pe. / a ^ ^^^^ 



36 ^ 12 \ho. 

A positive destabilizing term in (15T|) is derived from the first term in dHi"^^ /dz\z=i- 
We find that when /i is small, the trigger of the destabilization of the ice-water interface 
originates from the perturbed part of air temperature gradient at the water-air surface 
because —fifi^^\z=o of the second equation in fl33|) is proportional to fi when fi is small 
as shown in figure [3] (b). On the other hand, a negative stabilizing term in ( l5Ti) is 
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derived from the sum of second and third terms in dHi /dz\z=i. As shown in figure [3] 
(a), as /i increases, acquires non-zero values, hence the sum of second and third 
terms in dHi^^ /dz\z=i dominates over the first term and suppresses the instabihty. 
When fj^^^ly^^i ^ 0, we showed that there exists a phase shift of the water-air surface 
against the ice-water interface. This suggests that the instabihty and/or stabihty of 
disturbances of the ice-water interface is related to the magnitude of phase shift of the 
water-air surface, which will be discussed in the next section. As a result of competition 
between the first and second term in (1511) . we find from dai^^ /djj = that ai^^ acquires 
a maximum value at yU = [3{hQ/ay/ PeiY^^. From this, we obtain a simpler formula to 
determine the wavelength of ripples: A = 27r(a^/;.o-Pe//3)^/^ (for the dependence of A on 
Q/l, see Fig. 6 (a) in (Ueno 2007)). This formula includes two characteristic lengths a 
and ho- Indeed, using the typical values of a = 2.8 mm, Hq ~ 100 /im and Pei ~ 10, one 
centimeter scale wavelength is obtained from the above formula. It should be noted that 
this long-length scale is in contrast with the wavelength Ams = '^T^Vhdo obtained from 
the Mullins-Sekerka theory, which is of order microns (Mullins-Sekerka 1963). Here, Id 
is the thermal diffusion length and is usually a macroscopic length, whereas do is the 
capillary length associated with the solid-liquid interface tension, and is a microscopic 
length of order angstroms (Langer 1980, Caroli et al 1992). 

In the case of Ts\y=Q = T,|j^=^ = T^i + AT^i and Ti\y=^ = Ta\y=^ = Tia (Ueno 2003), 
the numerical results are shown by the solid lines in figures [5] (a) and (b). The dashed 
lines are the analytical results, (H7|) and ( l48l) . The deviation of the dashed line from 
the solid line in figure [5] (a) is mainly due to the neglect of the higher order of /i-Pe^. 
However, as shown in the inset in figure |5] (a), the analytical result is in good agreement 
with the numerical result as far as we are concerned with the long wavelength region 
of yU < 0.15. Here ^ = 0.15 corresponds to the wavelength of 3.8 mm for Q/l = 50 
[(ml/h)/cm] and 6 = 7i/2. 

In the case of Ts\y=(; = Ti\y=^ = Tsi and Ti\y=^ = Ta\y=^ = Tia + AT^^ (Ogawa 
and Furukawa 2002), the same ordinary differential equations were solved numerically 
with the same boundary conditions by replacing only the first equation in fl33|) with 
that in ( 1371) . and by neglecting the effect of the restoring force, i.e., a = in the 
last equation in (123|) . It is found in figure [5] (c) that there is a discrepancy between 
the analytical result (149|) (dashed-dotted line) and our numerical result (dashed line) 
calculated by us on the basis of their model but with no approximations. According 
to the stability analysis of the ice-water interface in the papers (Ogawa and Furukawa 
2002, Schewe and Riordon 2003), the instability of the ice- water interface occurs by 
the Laplace instability due to the thermal diffusion into the air. Moreover the flow 
in the thin water film makes the temperature distribution uniform, thus inhibiting the 
Laplace instability. They conclude that ripples of centimeter-scale wavelengths appear 
as a result of the competition between these two effects, and that the ripples on icicles 
should migrate downward. However, there are serious problems with this interpretation. 
First, our numerical calculation showed that even the length scale of ripples cannot be 
determined from their model because a* > for any wave number as shown by the 
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Figure 5. Non-dimensional amplification rate al and non-dimensional phase velocity 
Wp* versus non-dimensional wave number fi at Q// — 50 [(ml/li)/cm] and = 7r/2. (a) 
and (b) are in the case of ATgi ^ and ATia = 0. (c) and (d) are in the case of 
ATsi = and ATia ^ 0. Solid lines in (a) and (b): numerical results with boundary 
conditions in (Ueno 2003); dashed lines: analytical results (|47)) and (|48l) . Solid lines in 
(c) and (d) {a ^ 0) and dashed lines (a = 0): numerical results calculated by us with 
boundary conditions in (Ogawa and Furukawa 2002); dashed-dotted lines {a = 0): 
analytical results (|49| and ([50]). Here /i = 0.5 corresponds to the wavelength of 1.2 
mm. 



dashed line in figure |5] (c). This means that there exists no stabihzation mechanism of 
the ice-water interface. Second, according to the Laplace instability the latent heat is 
more rapidly lost from the convex surfaces than concave surfaces, resulting in faster ice 
growth on icicles's convex protrusions of the icicles than on its concave indentations. 
However, the Laplace instability cannot explain the translation mechanism of ripples. 

Our numerical results indicate that the approximation used to derive f j49|) is 
obviously incorrect. The same faulty approximation was made when deriving ( H9|) from 
our theoretical framework in (Ueno 2004). The comparison of the dashed line to dashed- 
dotted line in figure |5] (c) shows that the /i^ term becomes dominant for n > 0.13. 
Hence we cannot neglect the /x^ term in the equations and boundary conditions when 
deriving ( l49l) . In the case of (Ueno 2003), however, there exists already a stable region 
for fi < 0.13 as shown in figure |5] (a), and thus the long wavelength approximation 
neglecting the higher order of fi is valid. Even if we take into account the effect of a 
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in the model (Ogawa and Furukawa 2002), the situation is not improved as shown by 
the sohd hnes in figures [5] (c) and (d) as far as the boundary conditions in (1371) are 
used. The results are significantly different from figures [5] (a) and (b): the system is 
unstable for perturbations of any wave number and the sign of phase velocity is opposite. 
The leading cause of these differences in the two models originates from the boundary 
conditions between (133]) and fl371) when solving fl3T|) . 

4. Reconsideration of instability and stability of the ice-water interface 

From the mathematical expression indicated by the terms including a with minus sign 
in (HTll . it was suggested that the restoring force due to gravity and surface tension is 
an important factor for the stabilization of the ice-water interface on a long length scale 
of about 1 cm (Ueno 2003). However, the detail of the morphological instability and/or 
stability mechanism of the ice-water interface was not clarified. In the subsequent paper 
(Ueno 2004), it was shown that there exists a phase shift between a disturbed ice- water 
interface and the maximum point of heat fiux at its interface, and that the instability 
and/or stability of the interface is related to the magnitude of this phase shift. However, 
the cause of the occurrence of such a phase shift was not well understood. Here, this 
is investigated in detail by drawing the isotherm in the water film. In figures [6] (a), 
(b) and (f), the upper and lower solid hnes are the water-air surface and the ice- water 
interface, respectively. Figures [6] (c) and (d) show isotherms in the vicinity of the ice- 
water interface of figures [6] (a) and (b), respectively. 

It is convenient to express the temperature in the water film, T; = 7) + T/ = 
Tsi — Giy + HiGiC,^ in the non-dimensional form, as follows: 

Ti,{y.) = -J^^ = -y^ + 6My,)exp[at + ikx]. (52) 

-Lsl — -Lla 

We also express the temperature in the ice, Tg = Tg + T^' = Tgi + gs{y)exp[at + ikx], in 
the non-dimensional form. As far as kbo ^ 1 (6o ^ 1-6 mm for 1 cm ripple wavelength), 
can be approximated as follows: 

Ts^iy*) = J^"^ = 6bexp{fiy^){Hi\y^=o - l)exp[o-t + ikx]. (53) 

J- si — -L la 

At the ice- water interface y^, = C,*) taking the imaginary part for the perturbed part in 
(I52D and (ES]) yields 

TiAy,=c, = 7^.*ly*=C* = [{Hi'\y,=^ - 1)' + (i/f ^|,,=o)T^''^feW sin[fc(a; - V) + ©t.J, (54) 

where 9t^^ is a phase shift of the maximum point of temperature at y^. = against 
that of the ice-water interface (see the horizontal arrows in the upstream direction in 
figures E] (a), (b), (c) and (d)). On the other hand, at the water-air surface y^ = 
taking the imaginary part for the perturbed part in ( 152|) yields 

X 6b{t) sin[fc(x — Vpt) + OtjJ, (55) 
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Figure 6. (a)-(e) are in the case of ATsu 7^ and ATia* = 0. (f) is in the case of 
ATgi^ — and ATja, 7^ 0. (a) and (b) are isotherms in the water film at an unstable 
point (/X = 0.06, A = 9.6 mm) and a stable point (/i — 0.15, A ~ 3.8 mm) in figure [5] 
(a), respectively, for (5b = 0.05, Qjl = 50 [(ml/h)/cm] and Q — 7r/2. (c) and (d) are 
isotherms in the water film in the vicinity of the ice- water interface in (a) and (b), 
respectively, (e) 0^, (solid line): phase shift of the water-air surface; 0t^^ (dashed 
line): phase shift of the temperature at the ice- water interface; Ogj.-g^, (dashed-dotted 
line): phase shift of total heat fiux from the ice- water interface to the water and ice, 
against the ice-water interface. 9q,^_<j„ — ti/2 at = 0.092, which corresponds to 
the point a'"^^ = of the solid line in figure [5] (a), (f) is isotherm in the water film at 
fjL = 0.15 of the solid line in figured] (c). Qt^^ : phase shift of the temperature at the 
water-air surface against the ice-water interface. 



where B-r^^ is a phase shift of the maximum point of temperature at y^, = against 
that of the ice-water interface (see the horizontal arrow in the downstream direction in 
figure [6] (f)). 
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The conditions T^|y=c = ^i|y=C = '^si + AT^; in ([H]) and Ti\y=^ = Ta\y=^ = Tia in 
f lT^ adopted in (Ueno 2003) can be expressed in the non-dimensional form: 

Ts*\y,=Ct — ^/*|j/,=f, = AT^/*, TiJ\y^=^^ = Tat:\y^=(^^ = —1. (56) 

By comparing f l56|) to fl5^ and fl55|) . we obtain 

AT,,, = [{Ht^\y,=o - ly + {HP\y,=o)VW sin[fc(x - V) + ©t.J, 

\y*=i — -Jl Iy.=i5 ly.=i - \y*=i- [^') 

The last two equations in (!57|) are just the first equation in (!33|) . The dimensional 
form of the first equation in (lFr|) is AT^i = {Tgi — Tia)[{Hi^^\y^=Q — 1)^ + 
{Hi'^\y^=Q)^Y/^5bit)sm[k{x - V) + Using the solution of Hi at Q/Z = 50 

[(ml/h)/cm] and 6 = 7r/2, for which ai^^ acquires a maximum value at /i = 0.06, the 
maximum value of ATsi is 1.4 x 10~^ °C for the supercooling of Tgi — Tia = 0.03 °C of 
the water film (see Section [5]), fi = 0.06 and 5f,=0.05. The temperature deviation from 
Tsi = °C due to the Gibbs-Thomson effect evaluated at the same value of and Sb is 
of the order of 10~^ °C. Even if the value of AT,/ is extremely small but much greater 
than that due to the Gibbs Thomson effect, we cannot neglect the deviation because 
this contributes to the second terms in pij) . (142|) . ( 147|) and ( HHj) . which are represented 
by the dashed-dotted lines in figures S] (b) and (c). 

Figures |6] (a) and (b) show isotherms in the water film obtained from ( 152|) by 
using the solution Hi, determined by the boundary conditions in ( l33|) . It is found that 
the water-air surface is shifted by 6^, in the upstream direction against the ice-water 
interface, and that 9^^ increases as /i increases, as shown by the solid line in figure [6] 
(e). This phase shift is due to the effect of the restoring force in fi\y^=i in ( l39ll . The 
isotherm in the ice is determined by using f l53|) . Since the typical value of the thickness 
of the water film is about 100 /im, figures |6] (c) and (d) show isotherms around 10 /im 
from the ice- water interface. 

It should be noted that the isotherms in the water film are almost in phase with 
the shape of the water-air surface, as shown in figures |6] (a) and (b). Since the water- 
air surface is shifted in the upstream direction against the ice-water interface, the 
temperature distribution become non-uniform in the vicinity of the ice-water interface. 
This non-uniformity does not disappear even at the ice-water interface, as a result, the 
temperature at the ice- water interface deviates by AT^;*. Figures |6] (c) and (d) show 
that the maximum point of the temperature at the ice-water interface shifts by 0r^^ 
against that of the ice-water interface, which depends on n as shown by the dashed line 
in figure |6] (e) . For example, figure [6] (c) shows the isotherm near the ice- water interface 
with the wavelength of 9.6 mm. The deviation ATgi^, is positive on the upstream sides 
and is negative on the downstream sides of any protruded part of the ice-water interface. 
On the other hand, figure |6] (d) shows the isotherm near the ice-water interface with 
the wavelength of 3.8 mm. The deviation AT,;* is positive in any depressed region and 
is negative in any protruded region of the ice-water interface. 

We define the perturbed part of the non-dimensional heat flux from the ice- 
water interface to the water and from the ice to the ice-water interface, as g/* = 
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lm[-dTlJdy^\y,=(^^] and g,* = lm[-KfdTlJ dy^\y^=(^X respectively, where T/^ and T'^^ 
represent the perturbed terms in fl52|) and (l53l) . Hence, the total heat flux from the 
ice-water interface to the water and ice can be expressed as follows: 



qi* - Qs* = Sblia 



dif, 



y,=o 



+ Ki n{Hi\y^=Q — 1) f exp{at + ikx) 



(r) 



dy* 



y*=o 



+ 



dy* 



y*=o 



+ K^f,HP\y,=, 



1/2 



Sb{t) sm[k{x - Vpt) + eq^^_qj,{58) 



where 9„ 



is a phase shift of the maximum point of the total heat flux at = 



against that of the ice-water interface and this changes as fi increases, as shown by 
the dashed-dotted line in figure [6] (e). In order to avoid the temperature discontinuity 
by ATsi^ at the ice-water interface, the heat flux g^^, by the thermal diffusion occurs. 
The perturbed heat flux in the ice exists only in the vicinity of the ice-water interface 
because it is observed from ( 153|) that the non-uniformity of the temperature in the ice 
is exponentially attenuated far from the ice-water interface. 

When < 0g,^-g^» < vr/2 as is the case in figure |6] (c), the maximum point of 
qi* — Qs* shown by the vertical arrow is on the upstream side of any protruded part of 
the ice-water interface, which means that ice grows faster on the upstream side than 
on the downstream side of any protruded part of the ice- water interface. As a result, 
not only does the amplitude of perturbation grows, but ripples also move upward with 



time. On the other hand, when 7r/2 < 6q 



< TT as in the case in figure |6] (d), 



the maximum point of qu — Qs* shown by the vertical arrow is in any depressed region 
of the ice-water interface. This means that ice grows faster at any depressed part 
of the ice- water interface, and grows slower at any protruded part. Accordingly, the 
disturbance of the ice-water interface diminishes with time and such a disturbance 
eventually cannot be observed. We find that a phase shift between a disturbed ice- 
water interface and the maximum point of heat flux at its interface comes from the 
non-uniform temperature distribution at the ice-water interface due to the phase shift 
of the water-air surface against the ice-water interface. We also find that in order to 
explain the ripple migration it is necessary to cause an asymmetry in the temperature 
distribution between the upstream side and the downstream side of any protruded part 
of the ice- water interface. 

We define the characteristic time of shear rate as Tg^, which is just inverse of 
the shear rate 5*. The shear rate at the ice-water interface for the semi-parabolic 
shear flow Ui in f|T9|) is 5* = dUi/dy\y=o = 2uo//io- Hence, r^^ = l/(2uo//io) = 
[3{gsm6/h'iyQ/l]~^/^ is of the order of 10"^ s for the typical range of Q/l = 10 ~ 100 
[(ml/h)/cm] and 6 = 7r/2. We also define the thermal relaxation time ~ l/(/taA;^) and 
Ti ~ 1/(k;A;^) of fluctuations with wave number k, which are associated with the thermal 
diffusivity of the air, Ka = 1.87 x 10~^ m^/s, and that of water, ki = 1.3 x 10"''' m^/s. 



24 



respectively ( see (Langer 1980, Caroli et al 1992) for the relaxation time). For example, 
Ta and Ti are, respectively, of the order of 0.1 s and 10 s for about 1 cm ripple wavelength. 
It is convenient to introduce a characteristic wave number kc by Tsh = n. When Tsh < n, 
temperature fluctuations with a wave number smaller than kc are affected by the shear 
flow before they can dissipate thermally. On the other hand, when Tsh > ti, temperature 
fluctuations with a wave number greater than k^ can dissipate thermally without being 
affected by the shear flow (Onuki and Kawasaki 1979). In our system, the value of kc is 
about 10"^'^, which corresponds to a wavelength of about 200 /im for the typical values 
of Mo ~ 1 cm/s and ~ 100 /im. Therefore, the condition r^h ^ is satisfled at the 
ice-water interface with a wavelength of about 1 cm. 

Based on the two time scales mentioned above, we will explain why did we choose 
the boundary condition Ti\y=^ = Ta\y=^ = Tia at the water-air surface. We wfll also 
explain why the non-uniformity in the temperature distribution at the ice-water interface 
does not disappear, resulting in the temperature deviation ATgi from Tgi. A local 
temperature deviation of ATia from Tia at a disturbed water-air surface dissipates quickly 
by thermal diffusion in the air because shear stress is zero at the water-air surface. 
Therefore, it is reasonable to impose the boundary condition Ti\y=^ = Ta\y=^ = Tia 
at the water-air surface. Since Tgh Ta <^ ti for about 1 cm ripple wavelength, 
however, the temperature distribution in the water fllm is determined so as to adapt the 
instantaneous disturbed shape of the water-air surface satisfying the boundary condition 
Ti\y=£^ = Ta\y=^ = Tia before a local temperature deviation ATgi at the ice- water interface 
dissipates. This means that there is not enough time to relax the non-uniformity of the 
temperature at the ice- water interface thermally. That is why the temperature deviation 
ATsi from Tgi remains at the ice-water interface and the perturbed heat flux qs^, by the 
thermal diffusion is maintained in the vicinity of the ice-water interface in the ice. 

On the other hand, the conditions Ts\y=^ = Ti\y=(^ = Tgi in ffTTj) and Ti\y=^ = 
Ta\y=s^ = Tia + ATia in ( 1T5|) adopted in (Ogawa and Furukawa 2002) can be expressed 
in the non-dimensional form: 

Tg^:\y^=(^^ = Ti^,\y^=(^^ = 0, TJ^|j,_^=^^ = Ta*\y^=^^ = 1 ~l~ ATia*- (59) 

By comparing ( l59l) to ( 15^ and ( l55l) . we obtain 

n^'^h — 1 n^'-h — n 

-H; \y,=0 — i, -H/ \y,=0 — U, 

AT,,. = [{Ht\y,=, + fix=,r + {Hi%=, + fP\y,=irf'm sm[k{x - V) + e^.j. 

(60) 

The flrst two equations in fl60l) are just the flrst equation in fl37j) . Figure [6] (f) shows 
the isotherm in the water fllm obtained from ( l52i) by using the solution Hi, determined 
by the boundary conditions in (137|) . If we take into account the effect of the restoring 
force on the water-air surface in the model (Ogawa and Furukawa 2002), the water-air 
surface is shifted in the upstream direction against the ice- water interface as in flgure [6] 
(b). However, it should be noted that the isotherm in the water fllm is almost in phase 
with the shape of the ice-water interface, which is in contrast to the case in flgure [6] 
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(b). The maximum point of heat flux qi^ — indicated by the arrow is also different 
in between flgures[6] (b) and (f). In the case of flgure[6] (f), qs^ = because ATgi = 0. 
Under such a situation, if there exists a phase shift between the water-air surface and 
the ice-water interface, the non-uniformity of the temperature occurs at the water-air 
surface, as shown in figure [6] (f). The temperature at the water-air surface deviates by 
ATia from Tia- The maximum point of the temperature at the water-air surface shifts by 
Gj"^^ against that of the ice-water interface. Figure [6] (f) shows that ATia^. is negative on 
the upstream sides and ATia* is positive on the downstream sides of any protruded part 
of the water-air surface. The temperature distribution in the water film is determined 
so as to adapt the instantaneous disturbed shape of the ice-water interface satisfying 
the boundary condition Ts|y=^ = Ti|y=^ = Tgi before a local temperature deviation ATia 
at the water-air surface dissipates. However, this picture is physically inconsistent with 
the time scale Tsh -C ^ for the wavelength of interest. Actually, figure [5] (c) 
obtained from the boundary condition Ts|y=^ = Ti\y=^ = T^i shows that ripples with a 
characteristic length scale cannot be observed on the ice surface because all modes are 
unstable. 

The thermodynamics of fiuids under shear fiow is a challenging topic in modern 
non-equilibrium thermodynamics. For example, according to non-equilibrium molecular 
dynamics simulations of a system of spherical particles, the coexistence of crystal and 
shearing liquid flow cannot be accounted for by the equality of the chemical potentials 
of the crystal and liquid or by invoking a non-equilibrium analogue of the chemical 
potential (Butler and Harrowell 2002). In our system, the ice- water interface is in a 
non-equilibrium state under the influence of the boundary of the water-air surface as 
indicated in (!33l) and shearing water flow. Since such a non-equilibrium contribution 
would be expected to change the water chemical potential, there exists no physically 
reasonable deflnition of a non-equilibrium water chemical potential that would equal 
the chemical potential of the ice. Therefore, we did not impose the boundary condition 
Ts\y=(^ = Ti\y=(^ = Tgi at the ice-water interface, where T^i is the equilibrium freezing 
temperature only when the chemical potential of water equals that of ice. 

5. Experimental results 

In this section, theoretical predictions are compared with experimental results. As 
shown in flgure [TJ we used a pump which can control the water supply rate within the 
range of 50 to 500 ml/h. Water was pumped from the reservoir and dripped from the 
tip of the silicon tube at the top of a gutter on an inclined plane and of a round stick. 
A wooden plane with = 80 cm in length, / = 3 cm in width and 2 mm in thickness 
was inserted in the gutter of 2.5 cm in depth. The both sides and the back side of the 
gutter were covered with an insulation material to prevent the loss of latent heat at the 
sides. The stick was made of wood with dimensions of Ix = 80 cm in length and 6 mm in 
diameter. The thermal conductivity of wood is normally around 0.10 ~ 0.15 J/(mKs), 
which is much smaller than that of ice and aluminum. Since these instruments were set 
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Figure 7. Schematic view of apparatus set in a cold room. Water exits from the tip 
of the silicon tube covered with a band heater at the water supply rate of Q ml/h and 
flows down along a wooden plane and a stick. 



in a cold room, they were protected by a heating device in order to prevent the water 
from freezing in the sihcon tube. The temperature of the water dripping from the top 
at the rate Q ml/h was slightly above °C. The ceiling of the cold room was equipped 
with three fans. Although the fans were set to switch on and off periodically to make 
the temperature in the cold room uniform, large temperature fluctuations of ±3 °C 
around —9 °C were observed. The water reaches the supercooled state as it flows down 
along the plane and the stick. Ice grows from a portion of the supercooled water layer 
through which the latent heat of solidiflcation is released into the ambient air below 
°C. The rest of the water drips from the lower edge of the plane and the stick. Our 
measurement of the mean growth rate V of the ice produced on the 6-mm diameter 
round stick was 1.7 mm/h, which was almost independent of the water supply rate. 
This result is consistent with previous theoretical (Makkonen 1988) and experimental 
results (Maeno et al 1994). This is evident from (|30l) that V is independent of Q/l. 

The wavelengths in flgures |8] (a) and (b) are determined from the value of fi = kh^ 
at which ai^^ acquires a maximum value for a given Q/l and 6. Figure |8] (a) shows the 
dependence of the ripple wavelength on the angle of the inclined plane at Q/l =160/3 
[(ml/h)/cm]. As shown by the solid and dashed lines, our numerical and analytical 
results are in good agreement with the experimental results (A and • ). It is found 
that the wavelength of ripples increases with a decrease in angle. 

Figure IH](b) shows the dependence of the ripple wavelength on Q/l at 6 = tt/2. As 
shown by the solid and dashed lines, our numerical and analytical results show that the 
wavelength increases only gradually with an increase in Q/l. The experimental result 
(■) shows weaker dependence of the wavelength on Q/l than that expected from the 
numerical and analytical results, but the qualitative behavior and order of wavelength 
are almost the same. It should be noted that a portion of the supplied water freezes, 
and that the rest flows down the surface of the ice. Therefore, Q/l in Iiq, Pei and Rci 
should be replaced hy Q/l — {ps/ pi)Vlx from the mass conservation, where ps and pi are 
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sin e QjI [ (ml'li)/™] 



Figure 8. Wavelength of ripples of ice produced on a gutter, (a) The wavelength 
versus sm9 at Q/l =160/3 [(nil/h)/cm]. Solid line: numerical result; dashed line: 
analytical result (Ueno 2003); A: experimental result (Matsuda 1997); • : our 
experimental result, (b) The wavelength versus Q/l a.t 6 = tt/2. Solid line: numerical 
result; dashed line: analytical result (Ueno 2007); ■: our experimental result. 



the density of ice and water, respectively. Using the value of Ps/pi = 0.9, = 80 cm 
and V = 1.7 mm/h under the assumption that the ice is completely produced along the 
gutter from the top to bottom, unfrozen water is given by Q/l — 12 [(ml/h)/cm]. Hence, 
the values of h^, Pei and Rei are less than those estimated from Q/l supplied from the 
top. Therefore, the actual ripple wavelength is expected to be slightly less than what 
the solid and dashed lines show in figure [H] (b). 

Our theory also predict that the ripples move upward with |fp^,| = |fp|/V" ~ 0.6 at 
Q = 200 ml/h. This value is defined from jj, = kho at which ai^^ acquires a maximum 
value for a given Q/l and 6. It should be noted that the displacement of ripples depends 
on the growth rate V. Using the measured mean value of = 1.7 mm/h, we obtain 
\vp\ ~ 1 mm/h, meaning that the displacement over 4 hours is about 4 mm. Indeed, the 
observations in figures M (a) and (b) show that all ripples move upward. For example, 
ripples indicated by the arrows pass through the dashed lines. Although all ripples 
move upward, their speeds are not uniform because some ripples sometimes do not 
move when some portion of the ice surface is not covered with water. The measured 
mean displacements for 4 hours in figures [9] (a) and (b) are about 3.2 mm and 4.2 mm, 
respectively, which are of the same order as the theoretical results. In our experiment, 
Q from the top is kept constant. In the case of the ice produced on the round stick, 
the value of Q/l decreases as ice grows because the value of / increases with time t as 
27r(i?o + Vt) under the assumption that heat conduction into the wooden round stick 
through the ice is negligible, so that ice grows uniformly at V, where Ro is the stick 
radius. As a result, non- wetting parts on the ice surface increase as ice grows. The 
ripples produced initially almost disappeared over a 20-hour period due to sublimation. 
In this experiment, an upward movement of ripple was observed. This result is consistent 
with the observation that many tiny air bubbles are trapped in the upstream region of 
any protruded part of an icicle, and line up in the upward direction during icicle growth 



28 




Figure 9. A sequence of images showing upward movement of ripples of ice produced 
on (a) a 6-mm diameter round stick and (b) a gutter on a plane at 6 = 7r/2, after time 
6, 7, 8, 9, and 10 hours (from left to right) at Q — 200 ml/h. The mean displacement 
of ripples are (a) 3.2 mm and (b) 4.2 mm over 4 hours. 



(Maeno et al 1994, Ueno 2007). 

In the absence of heat conduction into the substrate through the ice, from the 
second equation in ( 12^ . we can estimate the degree of supercoohng of the water film at 
Tsi — Tia = LVho/Ki = 0.03 °C by using the measured value of = 1.7 mm/h and the 
typical value of ho = 100 fim. Infrared instrumentation was used to keep the surface 
temperature of a thin water film flowing on growing ice below °C (Karev et al 2007). 
The measurement showed that the surface temperature of the thin water film was always 
below °C. In our case too, it is necessary to measure the degree of supercooling of the 
water film accurately by such a non-destructive sensing technique. 

6. Summary and Discussion 

The validity of the approximations used in the two theoretical models (Ogawa and 
Furukawa 2002) and (Ueno 2003) was numerically investigated. There was an apparent 
discrepancy for the amplification rate between the analytical result in (Ogawa and 
Furukawa 2002) and our numerical result, in spite of solving the same governing 
equations with the same boundary conditions as those used in the model (Ogawa and 
Furukawa 2002). The characteristic length scale of ripples could not be determined 
under their boundary conditions. On the other hand, the author's analytical results 
(Ueno 2003, Ueno 2007) were in good agreement with our numerical results, and the 
theoretical predictions were confirmed in our own experiments: (i) the wavelength of 
ripples increases with a decrease in the angle of the inclined plane, (ii) the wavelength 
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increases only gradually with an increase in the water supply rate per width, and (iii) 
the ripples move upward. 

We also extended the theoretical framework (Ueno 2003) to include heat conduction 
into the substrate through the ice. If we take this into account, the wavelength depends 
on the ratio of the unperturbed temperature gradient in ice to that in water, , which 
includes the parameters such as ice thickness, thermal boundary layer thickness in the 
air, as well as substrate and ambient air temperature. All wavelengths of ripples on 
icicles observed in nature as well as those produced experimentally on a wooden plane 
and stick have a centimeter-scale. The thermal conductivities of these materials are 
small, and so the heat conduction into the substrate can be neglected. Hence, we can 
say that the kind of universahty of centimeter-scale ripples produced on an ice surface 
is due to Gl = 0. On the other hand, in the case of Gf ^ 0, our stability analysis 
predicts a critical value of Gf, below which the wavelength depends on the parameter 
Gf, but above which ripples do not appear. Our stability analysis would be applicable 
to prevent ripple formation on ice surfaces grown on a substrate for a given Q/l and 9. 

There is an analogy between wet growth and icicle growth (Makkonen and Lozowski 
2008). In the case of wet growth, water is collected from the impingement of supercooled 
water droplets. Whereas, in the case of icicle growth, water is supphed from melting 
snow and ice at the root of the icicle. In both cases, the surface of ice is covered with a 
supercooled water film, and ice grows from the portion of water film, by releasing latent 
heat into the ambient air below °C. From the point of view of engineering applications, 
the above analogy leads us to consider growth and morphology of ice on aircraft wings, 
wind turbine blades and aerial cables under wet icing condition, based on our theoretical 
framework. To do so, the present framework should be extended to include air flow, a 
supercooled water film motion driven by gravitational and aerodynamic forces, surface 
tension, and heat conduction through the ice into the object of cylindrical or arbitrary 
shape (Myers et al 2002a, Myers et al 2002b, Myers and Charpin 2004, Fu et al 2006). 
In our model, the flow of water fllm is driven by gravity only, and both surface tension 
and gravity act on the water-air surface of the water fllm flowing down an inclined plane. 
The present basic velocity proflle Ui^ in the water fllm was derived from the free stress 
condition at the water-air surface. If an aerodynamic force also acts on the water-air 
surface, the basic proflle of Uu would change from the half-parabolic form. How were 
our results modified by the effect of the aerodynamic force? The details will be discussed 
in a later paper. 
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